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1 Introduction 


Conventional spectral methods impose rigid requirements on the computational grids used. 
The grid points are the nodes of Gauss-like quadrature formulas (Gauss, Gauss Radau, or 
Gauss Lobatto (GL) formulas). These nodes are denser at the boundaries than in the middle 
of the domain. Although this property is suitable for boundary-layer problems, it may create 
difficulties for other types of problems, particularly those with disparate length scales that 
occur in multiple regions of the domain (e.g., diffusive burning or detonation and reacting 
mixing layers). The principle reason for the degradation in performance on these disparate 
problems is that the predetermined node points do not, in general, coincide with the features 
that are being resolved. Extensive mappings can concentrate the node points into regions 
more ideally suited for accurate resolution but present a serious limitation for complicated 
problems. For this reason, spectral multidomain techniques have an obvious advantage for 
complicated problems [1] — [3]. 

Another complication that conventional spectral methods have, is their implementation 
in complex geometries. Meshes that are predetermined present a significant constraint. 
Flexible mesh distributions are easily extended to geometries that are not tensor products 
of straight lines (to be shown in a later work). 

Spectral methods that are not constrained to specific nodal points would clearly be 
more flexible than conventional spectral methods. Specifically, a distribution of points that 
more closely approximates the disparate features in the domain could be adopted from 
the outset. Subsequent adaptation to solution features in the domain need not rely on 
smooth mappings. In addition, these “arbitrary-grid spectral techniques” could be used in 
conjunction with multidomain ideas. We focus on formalizing these ideas within the context 
of spectral techniques. 

In this paper, we present some ideas for constructing spectral methods with arbitrary 
grids. We demonstrate these ideas for a case of spectral solutions of hyperbolic equations; 
however, these ideas can be applied to any partial differential equation. To illustrate the 
basic idea, consider the following hyperbolic system of equations in conservation form: 


dU _ dFjJJ) 
dt dx 


-1 < x < 1 


( 1 ) 


with arbitrary initial and boundary conditions. For spectral methods, a polynomial (in the 
spatial variable x ) of degree N, U N (x,t), and a projection operator In are sought such that 


dl>N _ 0 In^{Un) _ o 

dt dx 


Of the spectral techniques, the most popular method is the Chebyshev collocation method, 
in which /^/(x) collocates f(x) at the Chebyshev GL points £,• = cos(^). Note that we 
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have here two projections; one involves the differentiation of F(Un), and the other involves 
the way that the equation is satisfied. Thus, the first application of the operator In occurs 
when we approximate qzP' by the derivative of the interpolation polynomial that interpo- 
lates F(U ) at the Chebyshev GL quadrature nodes. The second application of In occurs 
when we satisfy the approximate equation 


dU N 

dt 


dI N F{U N ) 

dx 


= 0 


at the Chebyshev GL points. 

The basic premise for unstructured spectral methods is that equation (2) does not have to 
be satisfied in the same manner in which the operation j § carr i e d out. In particular, 

the derivative operation can be carried out by interpolation at any selected points; the 
equation is satisfied by either a Galerkin formulation or by a collocation method at a different 
set of points. Most importantly, the equation must be satisfied correctly. 

Mathematically speaking, we can replace equation (2) with 


In 


3Un 

~dT 


8JnF(Un} 

dx 


= 0 


(3) 


where In 7^ Jn- 

In reference [4], a particular case with this approach has been discussed. The operator 
Jn was defined by the Chebyshev collocation operator, and In was the Legendre collocation 
operator. In the constant-coefficient case ( F(U ) = U), this method reduces to the Legendre 
collocation method with an efficient way of calculating the derivative by using Chebyshev 
collocation points. We now generalize this notion to an arbitrary set of points, which enables 
us to apply spectral methods in circumstances for which the grid points are not nodes of 
some Gauss quadrature formula. 

The method discussed in this paper is different from using a transformation to redis- 
tribute the grid points. The use of a transformation to redistribute the grid points involves 
approximation of the solution by a polynomial in the transformed variable; as a result, the 
approximation is not a polynomial in the original variable. Our method utilizes a polynomial 
in the original variable. Moreover, the new method can be applied to totally unstructured 
grids. 

Finally, it should be noted that the new method has many similarities with spectral 
elements, although the method of derivation is different. For instance, Patera [5] or Korczak 
et. al [1] used global polynomials (Lagrangian interpolants), passed through the Chebyshev 
collocation points, to obtain spectral elements. However, their work was not generalized to 
arbitrary grids. 
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2 The Differentiation Matrix for Unstructured Grids 


Consider the set of points (xo = 1, xi, X 2 , x^_i, xyv = —1), where the points xi, x 2 , xjv_i 
are arbitrary. Let /(x) be a function defined everywhere in [—1,1]. The interpolation 
polynomial //v(x) that collocates /(x) at the points Xj is given by 

f N {x) = j N f = J2 f( x i) L A x ) ( 4 ) 

i= o 

where the Lagrange polynomials Lj{x) are defined by 


L(x) 

L : (x) 

L 0 (x) 

Ln{x) 


(x - Xi)(x - X 2 )...(x - Xyv_ 2 )(x - XjV— 1 ) 

^ 1 < j < N - 1 

(1 — Xj)(x — Xj)L'(xj) - J - 

(1 + x)L(x) 

21(1) 

(1 — x)L(x) 

2X(— 1) 


(5) 

( 6 ) 

(7) 

( 8 ) 


The Lagrange polynomial evaluated at the discrete points xjt for k / j, is equal to 0; 
Lj(xk) = • 

We use as the approximation to Note that has two alternative repre- 

sentations; the first is obtained by differentiating (4) as 


d/jv(x) 

dx 


E f( x j) L 'j( x ) 


3=0 


(9) 


The second representation stems from the fact that dJN J r ^ is a polynomial of degree N — l- 
therefore, 


(io) 

dx j = o 

Equations (9) and (10) are used to relate the grid-point values of the derivative f'w{xj) 
to those of the function. The most obvious way is to equate the expressions in (9) and (10) 
at the grid points Xk (0 < k < N) to obtain: 

fN{ x k) = Y,L' ] (x k )f{x j ) (11) 

j=o 

To rewrite expression (11) in matrix form, we first denote 

/' = [/^(*o),...,/^(zv)] T , /= [/(x 0 ),...,/(xjv)] r 
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which yields 


( 12 ) 


/' = Df 

where the differentiation matrix D is given by 

D = (<*>>) = (13) 

Another method for expressing the equivalency between (9) and (10) is to state that the 
difference between these expressions (which is identically 0) is orthogonal to all polynomials 
of degree < N: 


ri JH 


I Y [f( x j) L 'j( x ) - f' N (xj)Lj(x)] L k (x)dx = 0 0<k<N 

J ~ 1 j = o 

The system of equations that follows from (14) can be rewritten as 

N N 

^ ] mkj f jy(*^', j) ^ ) ^k,j f i') 0 ^ k ^ N 

3=0 j=0 


where 


and 


m k ,j = J Lj(x)L k (x)dx 


s k,j = J i Lj(x)L k (x)dx 
In the matrix form, equation (15) becomes 

M /' = S / 

where 


M - (mi fcj) 0 < j, k < N 

and 


(14) 

(15) 

(16) 

(17) 

(18) 
(19) 


S = (a*j) 0<j,k<N (20) 

Equations (14) and (11) are different manifestations of the same fact: (9) and (10) are 
equivalent. Therefore, the differentiation matrices derived from (14) must be the same as 
the matrix derived from (11) (with the assumption that M is invertible): 


D = M -1 S 
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( 21 ) 



To prove this directly, we show that 


M D = S 


( 22 ) 


By writing the (i, k ) term on the left-hand side of (22), we obtain 

( MD )i,k = J2 m >A,k 

J-0 

If we substitute (13) and (16) into (2), then we get 

£ L :( x ) L 'k( x i) 

3=0 

We now use the fact that every polynomial of degree N is identical with its N -degree inter- 
polation polynomial. Thus, because L\.(x) is a polynomial of degree N — 1 and 

£ L'k(*,)L, (*) 

3=0 

is its interpolant at the points xj (0 < j < N) then 

Y.L i (x)L' t (x i ) = L’ k (x) 

3=0 

which yields 

(MB)* = Li(x)L' t (x) = s,, k 
(which is apparent from (17)). This establishes expression (22). 


(MB)* = /_' L.(x) 


□ 

Thus, we have defined a new method, based on the arbitrary distribution of points, to 
approximate the derivative of a function. The attractive features of the representation (21) 
of the differentiation matrix are summarized in lemma 3.1 and lemma 3.2: 


Lemma 3.1 : 

The matrix M defined in (16) is a symmetric positive-definite matrix. 
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Proof: 


The fact that M is symmetric follows immediately from the definition (16). In fact, 

m k j = J I Lj(x)Lk(x)dx = rrij t k 

We must show that M is positive definite. Let V be an N + 1 component vector: 


V = ( V 0 ,...,V N ) 


Then, 


V T MV = m i,3 V i V 3 
!=0 3=0 

Recall the definition of m,j from (16). We get 


_ _ ,i N N 

= / Y. VjLj(x) VjLj(x)dx > 0 

•'“l *•— n n 


(23) 


i=0 j = 0 

Clearly, the equality sign holds only if V is the null vector. 


□ 


Equation (23) can be interpreted in a different way. Let v(x) be the polynomial of degree 
N defined by 


v(ary) = Vj 0 < j < N 

so that 

N 

»(*) = Y1 V 3 L 3( X ) 
i= o 

Then, (23) can be rewritten as 

V t MV = J\ v{xfdx (24) 

Thus, every vector V can be identified with a polynomial u(:r) that takes the values of its 
components at the grid points Xj. The vector space norm 

v t mv 


is equivalent to the function space norm 



Next, we will consider the properties of the matrix S. 
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Lemma 3.2: 


Let S be defined in (17), and let V be defined as before. Then, 

V T SV = - «%) 


( 25 ) 


Proof: 


We start by showing that S is almost antisymmetric. ^From the definition (17) 

s k j = J L'-{x)Lk{x)dx 


and integration by parts, we get 


s k,j — Lj(l)Lk( 1) Lj( 1 )•£<*( 1) — sj t k 


We now use the definition of the Lagrange polynomials (6)-(8) and note that the bound- 
ary terms vanish for 0 / j, k ^ N to yield 


s k,j + = 4Ao — h,N^j,N 


Thus, 


N N 

V T SV = ££ VjVkSkj 

k=0j=0 

Y N N 

= o E T, VjV k (s kJ + S jlk ) 
z j= 0 

Y N N 

= VjVk(8kfl{>i,0 — 8k,N8j,N ) 

Z fc=0 j=0 



which completes the proof of (25). 


□ 
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As before, equation (25) has a natural interpretation in the polynomial space. Let u(:r) 
be the polynomial of degree N such that u(xj) = v r Then, 

N N 

v T sv = 

fc=oj=o 
N N 


= v i Vk I Lk(x)Lj(x)dx 

k =0 j =0 ■ /_1 

= j v(x)v ( x)dx 

= jb 2 (i) - » 2 (-i)] 


Note that 


v(l) = v 0 , u(-l) = v N 

Thus, (25) is an integration-by-parts formula. 

The last issue that we will discuss in this section is the relationship between differentiation 
matrices, based on different grid-point distributions. Consider two grids Xj and yj ( j = 
0 Let the Lagrange polynomial TJ(x) be defined as in (6)— (8), and let L^x) be 
defined in a similar way, based on the set of points yj. This defines two differentiation 
matrices (see (11)): 

d, = (<&) = [(if)'te)] 

and 

D, = (<&) = [(IJ)'(y;)] 

We now show that the two matrices are similar. 


(26) 

(27) 


Theorem 3.1: 


Define the matrix T by 


II 

ii 

'• 5 ' 

] (28) 

Then define 


(T-'k = [£»(*,)] 

(29) 

l 

D, = TD.T-' 

(30) 
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Proof: 


1. Because L y k is a polynomial of degree JV, 

E £'(*)££(*>) = i2(x) 

i=o 

If we substitute x = y m , then we get 

N 

X/ Lj {ym)L y k (xj) = Lj((y m ) = 
i=o 

which proves 1. 


2. Again, the Lagrange polynomials, based on the grid points y 0 , are polynomials of degree 
N, therefore, their derivative can be represented as 


WW*) = E £"(*;)(£*)'(*) 

i=o 

By the same token, 


(31) 


(L’)'(x) = Y.(L*)'(x,)Lf(x) 

/= 0 

Now, we substitute x = y m in (31) and (32) to get 

(£?)V) = EE i?fe)(£')'(*,)£f(fc) 

i=o /=o 


(32) 


(33) 


The left-hand side is the (m, i) element of D„ whereas the right-hand side is the (m i) 
element of T~ l D x T\ thus, (30) has been proved. 


□ 

3 The Legendre Galerkin Method Based on Arbi- 
trary Grids 

Consider now the linear form of (1): 

Ut{x,t) = U x (x,t) 

U(x,0) = f ( x ) 

U(l,t) = g(t) 
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- 1 < x < 1 


(34) 

(35) 

(36) 



We introduce a new method for the discretization of (34), based on the differentiation 
matrix introduced in the last section. Note that the differentiation matrix uses the arbitrary 
grid xj. With the new method, we seek a vector 

V, — [uo(f ), ..., )] 

that satisfies 

M^- - Su — Te 0 [u 0 — (37) 

at 

where 

e 0 = (1,0,0, ...Of 

The discussion on imposing the initial condition is deferred until later in the paper because 
of subtle issues that involve convergence. Here, we generally will not use 

Wj(0 ) = f{xj) 0 <j<N 

unless the grid points Xj have special properties. 

The structure of the matrices M and S, indicated in (34) and (37), leads immediately to 
the following stability result: 

Theorem 4.1 : 

The method described in (37) is stable for r > 

Proof: 

We multiply (37) by u T to get 

J,7 

u T M— - = u T Su — TU T e 0 [u 0 — </(t)] 
at 

(38) 

We use the symmetry property for M and the almost skew symmetric property (25) for 
S to obtain 

= ^(«o - u n) - TUo[u 0 - $(*)] (39) 

For stability, we consider the case g(t ) = 0; from this case we can clearly determine that 
if r > then 

u T Mu < 0 

2 dt 

and stability exists in the norm induced by the positive-definite matrix M. 
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( 40 ) 



□ 


The stability result (39) can be represented in a different way in view of the equivalency 
between vectors and polynomials established in (24). Specifically, let u N (x,t) be an JVth- 
degree polynomial such that 


u N (xj, t ) = Uj(t ) 0 < j < N 


Then, from (24) we see that 
1 d f 1 


L d [ t j \2 j 1 d _ T 

Wt L Mx ’ t) ix = 2dt u 


Mu 


= M -u 2 N )-rul 


— ^[^(M) 2 ~ u n{— l,i) 2 ] — TUw(l,t) 2 

Thus, for the polynomial u N (x,t) we have stability in the usual L 2 norm. 

Now, we examine equation (37) from yet another point of view. By, multiplying (37) by 
M _1 , we get 


du 


^ = M -1 Su — rM -I e 0 [«o ~ 9 {t)) 


or in view of (21), we obtain 


(41) 


du 

= Du-tM ^ofuo-^O]- 
The expression M 1 e*o can be evaluated explicitly. 

Theorem 4.2 : 

Let M be the mass matrix defined in (16). Define the residual vector r by 

M _1 e 0 = f = (r 0 , ...,r N ) T 

Then, 

r . ^N+ijzj) d~ 

3 2 

where Pn(x) is the Legendre polynomial of order N. 


(42) 


(43) 
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Proof: 


We must verify that if r satisfies (43), then the expression 

Mr = e*o 

is also satisfied. Substituting (16) into expression (3) yields 


N 


(Mr),- = 

3=0 


r 1 N 

= / L i {x)Y^L ] {x)r 1 dx 

J=0 

Substituting expression (43) into (44) yields 

(Mr),- = [' LiWjrLjix) 

J ~ l j=o 


(44) 


^N+lipj) ~b P N {*i) 
2 


dx 


(45) 


Because P^ +1 and P' N are polynomials of degree < N, they coincide with their iVth- 
degree interpolation polynomials; therefore, 

Y. L A x ) 

3=0 


p 'n+i( x j) + P 'n( x :) 


P v + i ( x ) + Pn( x ) 

2 


2 


so that 



Pn+i( x ) + M( x ) 

1 

-l 

2 


dx 


= Mi) 


' p N+1 (i) + p N (i) 


-Li(-l) 


Pn + i ( — 1) + Pn(~ 1) 


- J L\{x) 


Pn+ i(x) + Pn(x) 


dx 


Recall that 

Pn( l) = l,fiv(-l) = (-!)" 

and that P/v and pv +1 are orthogonal to all polynomials of degree < N; the last two terms 
in the right-hand side of (46) vanish, and we are left with 

(Mr"), = Li(l) = 6 it o 


which proves theorem 4.2. 


□ 
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Theorem 4.2 sheds a new light on the connection between the method defined in (37) 
that uses the arbitrary set of grid points xj and the Legendre Galerkin method. They are 
the same method. 

Theorem 4.3: 


The method defined in (37) is equivalent to the Legendre Galerkin method. 


Proof: 


Define 


N 


j = 0 


where Uj(t) are the elements of u defined in (37). Then, v,N(x,t) satisfies the error equation 


du]v(x,t ) du^(x,t) 


dt 


dx 


— T 


P > N+li x ) + Pn( x ) 


[Ml, t)-g(t)} 


(46) 


The error equation is satisfied because both sides of expression (46) are polynomials of degree 
N. The two sides agree at A T + 1 points xj ( j = 0 by virtue of (37), which indicates 
that they are equivalent. Because the right-hand side is orthogonal to all polynomials of 
degree N that vanish on the boundary x = 1, this error equation is the same equation that 
is satisfied by the Legendre Galerkin method [6]. 


□ 


As equation (46) demonstrates, the precise method for imposing the boundary conditions 
affects the overall behavior of the method. Section 3 shows that two differentiation operators 
defined on different grids are similar and, thus, have the same eigenvalues. We now show 
that the modified differentiation matrix also has this property. Equation (42) produces a 
modified differentiation matrix (i.e., a differentiation matrix that takes into account the 
boundary conditions): 

D-tR 

where the boundary matrix R is defined as 

R\ ,j — (4 7) 


13 



Suppose now that we have two grids Xj, y 3 (j = 0, N). We have shown in theorem 2.3 
that D x and D y are similar: 

D y = TD x T~ l 

where the matrices T and T -1 are defined in (28) and (29). We show now that the same 
similarity transformation exists for the modified differentiation matrices. That is, 


D y - rR y = T ( D x — tR x ) T 


-l 


or (with theorem 3.1) 


R y = TR x T~ l 

Consider element (i, j) of the right-hand side: 


(TR t T-'). = ££T, 

,,J i=0 m= 0 

= EE Lf( yi )r,6o„L)M 

1=0 m = 0 


(48) 


We recall that 


r i = 


^AT +1 ( a: 0 + P'n( x i) 
2 


= Rn{xi) 


where Rn ( x ) is a polynomial of degree N and is, therefore, equal to its interpolant. Thus, 


(TR x T~ l ) t . = R n ( yi)L V j{x 0 ) = R N {yi)6 jft 

which proves that the similarity transformation is valid even for the modified derivative 
matrix. 

The Legendre Galerkin method defined by equation (37) is stable; therefore, the initial 
error is not amplified. However, the effects of initial conditions must be carefully taken into 
account. We know that polynomials based on arbitrary grid distributions may generally be 
nonconvergent (the Runge phenomenon). 

The initial error can be decreased with the number of mesh points N by constructing the 
Chebyshev interpolation as an initial condition. Thus, let 


ii = cos(^) 0<j<N 

(49) 

L'(z) = (1 - x 2 )T^(x) 

(50) 

L ( (x) _ 

(51) 


The Chebyshev approximation for the initial condition is, then, 

ttv/M = £/«,)!$(*) 

J=0 
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so that the recommended initial approximation will be 

/(*,) ~ £/(&)£<(*/) 

3=0 

This approximation will provide a convergent approximation for the initial condition. Of 
course, the Chebyshev approximation is not the only possibility; any other spectral or pseu- 
dospectral approximation would do as well. 

We now briefly discuss the issue of implementation. Two methods are available for 
implementing the arbitrary-grid spectral methods. The first method is to form the matrices 
M and S by carrying out explicitly the integrations in (16) and (17). (This technique is 
utilized in the two examples presented later in the text.) This procedure is done once and 
for all for every given set of grid points. Then, the equations are solved as described in 
(37). A more convenient method that does not involve evaluating integrals is to use the 
differentiation matrix D defined in (13) and solve the system (42) with the identity 

jyj-i ^ _ ^jv+iC^i) T ^/v( 3 'j) 

proven in theorem 4.2. For a large N , the method that will be the most successful is the one 
with the least sensitivity to round-off errors. This point has not been fully investigated at 
this time. 

Finally, an observation in regard to the maximum allowable time step for the arbitrary- 
grid spectral schemes. All spatial operators have the same eigenvalues, regardless of the 
spatial distribution of points (48). Therefore, the maximum allowable time step is the 
same for all schemes. Stability is a matrix property, and depends on all the points in 
the distribution. This observation is somewhat counter to the conventional finite-difference 
notion, in which the maximum time step is governed by the smallest grid spacing. 


4 The Legendre Collocation for Unstructured Grids 

The Legendre collocation for unstructured grids involves the approximation of the integrals 
in (16) and (17) by the Gauss-Lobatto-Legendre (GLL) quadrature formula. Let ( r) 0 = 
liVu-'tilN-uVN = -1) be the nodes of the GLL quadrature formula and u>/, 0 < / < N be 
the weights. We define a new mass matrix M c by 

N 

= £ Lj{T)i)L k (vi)u3i (53) 

1=0 

where the Lj{x) are the Lagrange polynomials at the points ( x 0 = l,x 1 ,x 2 ,...,x N ^ 1 ,x N = 
— 1). Note that this is an arbitrary set of grid points. 
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The matrix M c may be different from M because the GLL formula is exact to order 2 N— 1 
and Lj(x)L k (x ) is a polynomial of order 2 N. The matrix M c is, however, a symmetric and 
positive-definite matrix. 

By introducing quadrature to equation (17), we define a new stiffness matrix S c as 

s c(i,j) = ^2 L 'j(vi)Lk(m) u i ( 54 ) 

1=0 

Note that because of the exactness of the GLL formula, the sum on the right-hand side of 
(54) is the same as the integral in the right-hand side of (17); therefore, 

S c = S 

For this reason, the property (25) is true for the stiffness matrix S c also. 

The uniqueness of the differentiation matrix D also yields 

M?S e = M _1 S 

which does not contradict the fact that 

M C ^M 

because the matrices S c and S are singular. 

In the Legendre collocation method of (34) with arbitrary grids, we seek a vector 

u = [u 0 {t),...,u N (t)] T 

that satisfies 

clii 

M c — = S c u - re 0 [u 0 - g{t)] (55) 

where 

e 0 = (1,0,0,. ..Of 

Alternatively, 

fiv 

^ = Du-rM^eoWo-git)] (56) 

The stability of (55) follows immediately from the fact that M c is symmetric positive 
definite and S c satisfies (25). Our aim is to show that (55) is equivalent to the usual Legendre 
collocation method. 
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Theorem 4 . 1 : 


Let M c be the mass matrix defined in (53). We define the residual vector f by 

M c - X e 0 = r = (r 0 ,...,r N ) T 

Then, 

rj = P N ( Xj )(l + Xj )^ (57) 

where Pn(x) is the Legendre polynomial of order N. 

Proof : 

We start by noting that the nodes „ of the GLL formula are the zeroes of the polynomial 

P'n(x)( 1 - 1 1 ) 

Also, because P^(l + x) is a polynomial of degree N , 

N 

E L : (x)P' n {x } ){ i 4- Xj ) = P' N (x){ 1 + x) 

J= o 

Therefore, 

2N 7 = j:M c (i,j)P' N ( Xj )( l+Xj ) 

j—0 
N N 

= E E L j(Vi)Li(r),)P' N ( Xj )( 1 -(- * -) W/ 

j=o 1=0 

N 

= E ^( t u)Pn( t Ii)(^ + 

/=0 

= 2 jV 2 «V t - 

which proves the theorem. 



( 58 ) 


such that 


duyjxkJl = ( 59 ) 

dt j-Q 

where 

/? n (x) = (1 + x)F n (x) 

This approach is equivalent to the Legendre Collocation Method [6]. 

The extension of the arbitrary-grid Legendre collocation method from the linear case (34) 
to the solution of the nonlinear case (1) is immediate. The issue of implementation could be 
significant. To avoid computing the points r, h the best choice is to use the formulation (56) 
rather than (55). In this case, M c and S c do not need to be computed. 

At this stage, note that for the case 

Rn{x) = Pn ( x ) 

we have the Legendre Tau method, with the additional property of an improved time step. 
However, we do not have the representation of the Legendre Tau method in the form of (55). 


5 Unstructured Grids for Unbounded Domains: La- 
guerre Methods 


Consider the equation 


6U_ _ _dU_ 

dt dx 

U{0,t) = g{t ) 

U{x, 0) = h(x ) 


0 < X < oo 


Note that the domain is semibounded. Note also that if g{t) - 0, then 
— r e~ T U 2 {x,t)dx = -f e~ T U 2 (x,t)dx 

dt Jo 

Assume that we have an arbitrary set of grid points 


(x 0 = 0,Xi,...,Xjv) 


(60) 

(61) 


(62) 


In the Galerkin procedure, we approximate the derivative of a function f{x) whose values 
at Xj are given by the derivative of its interpolant f N {x). After we define 

L(x) = (x- x 0 ).-{x - x N ) 
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( 63 ) 


(64) 


we define the Lagrange polynomials by 

Lj(x) = 


L(x) 


(x Xj)L ( Xj ) 

The derivative of the interpolant f N (x) has two equivalent expressions: 


dMjz) 

dx 


N 


= E/(*i) 


and 


i=o 


N 


dLj(x) 

dx 


(65) 


d M x ) _ r/, u , . 

^ J N\ x j)Lj{x) 

j= o 


In the Galerkm Laguerre method, we express the equivalency between the 


( 66 ) 


expressions by 


J yoo N 

J = 0 


f{Xi) dT - foMW*) 


L k (x)dx = 0 0 < k < N (67) 


Equation (67) defines the differentiation matrix D. In fact, if we define 


= {Lj,L k ) 


( 68 ) 


and 


Sk,j — (Lj,Lk) 

where the scalar product (u, v) is defined 


(69) 


as 


then we get 


roc 

(u, v ) = e~ x u(x)v(x)dx 


D = M _1 S 


As before, the differentiation matrix is unique. The manner in which the matrices M and S 
are constructed leads immediately to the following lemma. 


Lemma 6.1: 


The matrix M is symmetric positive definite. The matrix S satisfies 

S + S r = M — diagonal(l,0, 0, ...,0) 


(70) 
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Proof: 


iFrom the definition of the matrix S, we have 

stj = (L'jM 


fOO 

= / e~ x LAx)Lk{x)dx 

JO 

roc , 

= -Lj{0)L k (0) - J o e~ x Lj(x)L k (x)dx 

fOO 

+ / e~ x Lj(x)Lk(x)dx 
Jo 

By using the definition of the matrix M and the properties of the Lagrange polynomials, we 
get 

Sfcj = — — s j,k + m k,j 


which proves (70). 


□ 


To discretize (60), we introduce the unknown vector 

U = [tf 0 (*)>— , u N{t)} T 

that satisfies 

m- = -sa-T«y«„-s(i)] < 72 > 

dt 

The stability is immediate, as shown in the following lemma. 

Lemma 6.2 : 

Let u satisfy (72), with g(t) = 0. Then, we have the energy estimate 

-u T Uu = -u T Uu - {2 t - l)u 2 0 (73) 

dt 


Proof : 

Equation (73) follows immediately from multiplying (72) by u T and using (70). 

□ 
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Lemma 6.1 implies that the method is stable, provided that r > Note that the energy 
estimate (73) for the approximation is nearly the same as for the differential equation (64). 

We still must show that the method described in (72) is equivalent to the Laguerre 
Galerkin method. We begin by rewriting (72) as 


~j[ ~ ‘Su-rM 1 e 0 [u 0 -g{t)\ 

The key issue is to identify the vector 

M _ 1 e 0 

which is done in the following theorem. 

Theorem 6.1 : 

Let M be the mass matrix defined in (68). Define the residual vector r by 

M - 1 e 0 = f = (r 0 ,...,r N ) T 

Then 

r ■ - jLr(°) I 

3 ~ dx N+1 '*=*> 

where £$ is the Laguerre polynomial of order N. 


(74) 


(75) 


Proof: 


We must verify that f satisfies (75) and that 

Mr = e 0 

We begin by expanding (Mr) as 


N 


(Mf) t . = 


3 = 0 

TOO 


N 


j e x Li{x)Y^Lj{ x ) r jdx 

0 7=o 


( 76 ) 


If we substitute (75) into (76 ), then we get 


Because j- x C §\ , is a polynomial of order N, it coincides with its interpolant; therefore, 


EM*) 


J=0 



( 0 ) | 
N+l 


-1 

dx 



(x) 


Thu s , ^ , 

(Mr) ; = jf°° /^(xjdx 

If we integrate the right-hand side by parts, we get 

(Mr), = -L.XOJCgirtO) + jT «‘'£i( I ) £ S?+i(*)* ! 

The last two terms on the right vanish because of the orthogonality of £ (0) , and the first 
term vanishes if i ^ 0; thus, 

(Mf)i = -ho 


and the theorem is proven. 


□ 


Another method for getting the Lageurre method on the grid x 2 is to seek a polynomial 
«jv(x,<) such that 


o,o- 1/(0] ( 78 ) 

dt j =0 

where ( LG) N (x ) is the A r th- degree Laguerre polynomial. This approach is the Laguerre 
collocation method. 


6 Numerical Results 


We now test the previous theoretical results with two numerical examples. The linear equa- 
tions (34)— (36) are solved with f(x) = sin(Trx), g(t) = sm[ir(l + f)]> and the exact solution 
U(x,t) = sin[ 7 r(x + t)}. A variety of grids, from Chebyshev to “randomly generated grids, 
are used to test the accuracy and stability of the method. For all calculations, 128-bit 
arithmetic is used to ensure adequate precision. 
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Figure 1 shows the refinement study on five different grids: 


1. Uniform grid x 3 = (y _ 0 , N) 

2. Chebyshev grid x, = cos(^) 

3. A linear combination of the uniform grid and Chebyshev grid (i e x = 0 5 2j ~ n 4. 

0.5cos(f)) ' " ^ 

4. Chebyshev 2 (i.e., x 3 = cos 2 (|/)) 

5. (Chebyshev)’ for -1 < * <0 and (Chebyshev)i for 0 < * <1, where 

(Chebyshev)* is defined by the grid points Xj = cosi(^). 


The log 10 of the £ 2 error, plotted against the number of points in the approximating 
polynomials is shown in Figure 1 . The problems are run to the physical time T = 2. The 
convergence is exponential for all cases until machine round off is encountered. These results 

are consistent with the previous numerical results. (Note that the Chebyshev grid is the 
least sensitive to round off.) ° 

The Legendre Galerkin method defined by equation (37) is stable; therefore, the initial 

error is not amplified. However, the effects of initial conditions must be carefully taken 

into account. We know that polynomials based on arbitrary grid distributions generally 

may be nonconvergent. This property, called the Runge phenomena, is easily demonstrated 

by approximating the function f(x) = rffer (-1 < * < 1) on a uniform grid. The 

global approximating polynomials oscillate wildly at each end of the domain, which yields 

a poor approximation in those regions. The Runge phenomena is alleviated by using a grid 

distribution (like the Chebyshev grid distribution), which clusters points near the boundaries 
—1 and 1. 

Figure 2 illustrates that a Runge-like phenomena exists within the arbitrary-grid spectral 
methods if special precautions are not taken in the initialisation step. In this problem, the 
linear equations (34)-(36) are solved with ffx) = — J. a (f) _ x , 

solution U(x t) — 1 cp, . . / i+W*)P’ ~ i+RT+ijF’ and the exact 

solution u(x,t) - T+iTc^+oi* ' The simulation is run to time T = 0.001 (a physical time that 

occurs well before the influence of the initialization is lost.) (Running to a physical time 

> 2 yields exponential convergence on all grids.) Convergence is achieved only for the 
Chebyshev grid distribution. 


The source of the error in this problem is the failure of the arbitrary grid that approx- 

unates the polynomial to converge to the initial condition. For small times (less than 1 

convective sweep), erroneous information is left in the domain, and the resulting method is 

nonconvergent. By changing the problem slightly, however, convergence can be recovered on 
all grids. 
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To initialize the problem, we must construct an approximation to the mitral condition 
/(x) based on the grid points x, (0 < j < N). We want to keep the flex.b.hty and ng. 

structure of the original grid distribution; however, the interpolat.on polyno^al 

the grid points generally is not convergent. Therefore, we use the method 
(49) and (52). With this initialization, spectral convergence is recovere . 

7 Conclusions 

A new technique for implementing spectral methods for hyperbolic equations has been de«l- 

this reason, this method is referred to as an arbitrary-gnd spectral method. Both Galerk 
and collocation formulations are presented for the conventional Legendre method, and 
Galerkin formulation is presented for the conventional Laguerre method 

The basis for the stability of the unstructured spectral schemes relies - ^ 
energy norm in all cases. Stability is proven for the constant coefficient hyperbolic sy - 
All unstructured spectral methods utilize a “weak” imposition of the boundary con 
similar to the technique used in the penalty formulations of the finite element method, 
this imposition, the complete differentiation matrix, including boundary cond, turns, » sim.lar 
to (i.e., it has the same eigenvalues) the conventional different.at.on operator; therefore, 

matrix behaves similarly. arViitrarv- 

The new formulations are demonstrated on two scalar hyperbolic problems. The arbitrary 

grid Legendre Galerkin method is used in both cases. Exponential accuracy is shown in both 

cases on arbitrary grids. Care must be exercised in the initialization proce ure o ensure 

convergence of the new schemes. 
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